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comments. 


Consider  linear  programming  problems  of  the  form 


Minimize  {c,x)  (P) 

subject  to  Ax  =  b,  X  >  0, 

where  c  is  an  m-vector,  A  is  an  nxm  matrix,  b  is  an  n-vector,  and  (•,  )  denotes  the  usual  Euclidean 
inner  product.  In  our  notation,  all  vectors  are  column  vectors  and  superscript  T  denotes  the 

transpose.  We  will  denote  by  91"'  (91")  the  m-dimensional  (n-dimensional)  Euclidean  space. 

For  any  vector  x  in  9?"*,  we  will  denote  by  xj  the  jth  component  of  x.  For  any  positive  vector  x 

in  91"’,  we  will  denote  by  the  mxm  positive  diagonal  matrix  whose  jth  diagonal  entry  is  the  jth 
component  of  x.  Let  X  denote  the  relative  interior  of  the  feasible  set  for  (P),  i.e. 

X  =  {  x€9("’  I  Ax  =  b,  X  >  0  }. 

We  will  also  denote  by  e  the  vector  in  91"’  all  of  whose  components  are  I's.  "Log"  will  denote  the 

natural  log  and  IMIj,  11112  will  denote,  respectively  the  Li-norm  and  the  L2-norm.  We  make  the 
following  standing  assumption  about  (P): 

Assumption  A: 

(a)  Both  X  and  {  ue  91"  I  A’^u  <  c  )  are  nonempty. 

(b)  A  has  full  row  rank. 

Assumption  A  (b)  is  made  only  to  simplify  the  analysis  and  can  be  removed  without  affecting  either 
the  algorithm  or  the  convergence  results.  Note  that  Assumption  A  (a)  implies  (cf.  [3],  Corollary 

29.1.5)  that  the  set  of  optimal  solutions  for  (P)  is  nonempty  and  bounded.  For  any  e  >  0,  consider 
the  following  approximation  of  (P): 

Minimize  f£(x)  (P^) 

subject  to  Ax  =  b, 


2 


where  we  define  f£:(0,‘>o)"'-49?  to  be  the  penalized  function: 


W  =  <c,x>  -  eXj  log(Xj).  (1.1) 

Note  that 


l/xj 

l/(Xj)2  ...  0 

Vfg(x)  =  c  -  e 

• 

V2f^(x)  =  e 

: 

1/Xjj, 

h.  ^ 

0  ...  l/(x^)2 

to  M 

The  literature  on  Karmarkar's  algorithm  [1]  is  too  vast  to  survey  and  we  will  not  do  so  here.  Our 
approach  to  solving  (P),  which  is  similar  to  that  taken  in  l5]-[6],  [1 1],  is  to  solve,  approximately,  a 

sequence  of  problems  {(Pjr)),  where  {e*”)  is  a  sequence  of  geometrically  decreasing  scalars.  The 
approximate  solution  of  (Pgr),  denoted  by  x*",  is  obtained  by  solving  the  quadratic  approximation  of 

(Pgr)  around  x^*'.  The  novel  feature  of  our  algorithm  is  its  simplicity,  both  in  the  description  and  in 
the  analysis.  And  yet  it  has  an  excellent  complexity.  Also,  our  algorithm  is  unusual  in  that  it  is  a 
primal  affine-scaling  algorithm  that  generates  feasible  primal  and  dual  solutions. 


This  note  proceeds  as  follows:  in  §2  we  show  that,  given  an  approximately-optimal  primal  dual 

pair  of  (Pg),  an  approximately-optimal  primal  dual  pair  of  (Pot),  for  some  ae  (0,1),  can  be  obtained 
by  solving  a  quadratic  approximation  of  (Pg).  In  §3  and  §4  we  present  our  algorithm  and  analyze  its 

convergence.  In  §5  we  discuss  the  initialization  of  our  algorithm.  Finally,  in  §6  we  give  our 
conclusion  and  discuss  extensions. 


1  Accasslon  For  y  \ 
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and  let  x  be  any  element  of  X  and  let  u  be  any  element  of  91".  We  replace  the  objective  function  by 
its  quadratic  approximation  around  x  =  x.  This  gives  (cf.  (1.2)) 

Minimize  (c-€(D)'le,z)  +  e(z,(D)‘2z)/2 
subject  to  Az  =  0, 

where  we  denote  D  =  Dj.  The  Karush-Kuhn-Tucker  point  for  this  problem,  say  (z,u),  satisfies 


c  -  e(D)‘ie  +  e(D)‘2z  -  A'l'u  =  0, 
Az  =  0. 


(2.1a) 

(2.1b) 


Letd  =  (D)'lz.  Solving  for  d  gives 


ed  =  [I  -  (AD)T(AD2AT)-lAD]r, 


(2.2) 


where  we  denote 

7  =  -Dc  +  ee  +  (AD)'nj. 

Note  that,  since  the  orthogonal  projection  is  a  nonexpansive  mapping  (with  respect  to  the  I^-norm), 
we  have  from  (2,2) 


Ildil2  ^  IF  112/e. 


(2.3) 


Let 


X  =  x+z. 


(2.4) 


D  =  D^,  and  A  =  Dj.  Then  D  =  D  +  A^  and  hence 


-Dc  +  ee  +  (AD)’l'u  =  [-Dc  +  ee  +  (AD)Tu]  +  A[-Dc  +  (AD)’^u] 
=  ed  +  A[-Dc  +  (AD)Tu] 

=  A[ee-Dc  +  (A^)Tu] 
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=  eAd, 

where  the  second  and  the  fourth  equality  follow  from  (2.1a).  This  implies  that  (with  dj  denoting  the 
jth  component  of  d) 

ll-Dc  +  ee  +  (AD)Tull2  =  ellAdll2 

<  ellAdlli 

=  eXj(dj)2 

=  e(lldll2)2 

<  (irFll2)2/e,  (2.5) 

where  the  first  inequality  follows  from  properties  of  the  Lj-norm  and  the  I^-norm  and  the  second 
inequality  follows  from  (2.3). 

Consider  any  jJe  (0,1)  and  any  scalar  a  satisfying 

(|32+ml/2)/((J+ml/2)  <  a  <  1.  (2.6) 

Let  e'  =  cx£  and  r'  =  -Dc  +  e'e  +  (AD)’*'u.  Then 

Ilr'll2/e'  =  ll-Dc +  txee+(AD)T'ull2/((xe) 

<  ll-Dc  +  ee  +  (AD)’*'ull2/((xe)  +  (l-a)  m*/2/a 

<  (lFll2/e)2/a  + (l/a-l)  mi/2^ 

where  the  first  inequality  follows  from  the  triangle  inequality  and  the  second  inequality  follows  from 
(2.5).  Hence,  by  (2.6),  if  Ifrll2/e  <  p,  then  Ilr'll2/e'  ^  p.  Funhermore,  by  (2.3),  we  have  that  Ildll2  < 
P  <  1.  Hence  e+d  >  0  and  (cf.  (2.4))  x  >  0.  Also,  by  (2. lb)  and  (2.4),  Ax  =  A(x+z)  =  b. 

For  any  e  >  0,  let  pg:(0,<»)"'x91"— >[0,'»)  denote  the  function 
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Pe(y.p)  =  ll-DyC  +  ee  +  (ADy)T'pll2/e.  (2.7) 

We  have  then  just  proved  the  following  important  lemma: 

Lemma  1  For  any  E  >  0,  any  pe  (0,1 )  and  any  (x,u)e  Xx5l"  such  that  p£(x^)  <  3,  we  have 
(x,u)eXx9?",  Pa£(x,u)  <  p, 

where  a  =  (pVm’^)/(P+m^/2)  and  (x.u)  is  defined  as  in  (2.1a),  (2.  lb),  (2.4). 


3.  The  Homotopy  Algorithm 

Choose  a  =  (l/4+ml/2)/(i/2+ml/2).  Lemma  1  and  (2.1a)-(2.1b),  (2.4)  motivate  the  following 
variant  of  Karmarkar’s  algorithm,  parameterized  by  two  positive  scalars  y  <  E: 

Homotopy  Algorithm 

Step  0:  Choose  any  (x^u')e  Xx5?"  such  that  p£(x^u0  ^  1/2.  Let  e’  =  E. 

Step  r:  Compute  (z‘’+l,u‘‘+')  to  be  a  solution  of 


Er(Dxr)-2  -at 

z 

E^OxO'e  -  c 

A  0 

u 

1-  • 

0 

Set  x''+l  =  x*^  +  z'’'*'^,  =  OE*'. 

If  E*^^^  <  y,  terminate. 

[Note:  For  simplicity  we  have  fixed  P  =  1/2.]  We  gave  the  above  algorithm  the  name  "homotopy" 
because  it  solves  (approximately)  a  sequence  of  problems  {(Pgr))  that  approaches  (P)  (see  [2]). 


By  Lemma  1,  the  homotopy  algorithm  generates,  in  at  most  (log(Y)-log(e))/log(a)  steps,  an 
(x,u)e  satisfying 

ll-D^c  +  “ye  +  (ADx)Tull2  ^y/2,  (4.1) 

Ax  =  b. 

Since  Y>  0,  (4.1)  implies  that 

0  <  D^c  -  (ADx)Tu  <  (3Y/2)  e. 

Since  is  a  positive  diagonal  matrix,  multiplying  both  sides  by  (D^^)  *  gives 
0  <  c  -  ATu  ^  (3Y/2)  (Dx)-ie. 

This  in  turn  implies  that,  for  each  j€  { l,...,m}. 


0  <  Cj  -  (Aj,u)  <  3y''^/2  if  Xj  >  y'^, 

0  <  Cj  -  (Aj,u)  otherwise, 

where  Cj  denotes  the  jth  component  of  c  and  Aj  denotes  the  jth  column  of  A  (note  that  an  analogous 
argument  shows  that  u*"  is  dual  feasible  for  all  r).  Also,  since 

log(l-6)  =  -5  -  §2/2  -  §3/3  -  §'‘/4  -  . . . 

^  -§(1  +  (§/2)  +  (§/2)2  +  (§/2)3  +  . . .) 

=  -§/(l-§/2)  =  -l/(l/§-l/2), 

for  any  §€  (0,1),  we  have  that 


log(a)  =  log(l-(244m*'2)-i) 


^  -l/(2+4m»/2-i/2). 


Hence  we  have  just  proved  following: 

Lemma  2  For  any  positive  scalars  y^t,  the  homotopy  algorithm  generates,  in  at  most 
(3/2+4m*/2).(iog(e)-log(Y))  steps,  a  pair  of  optimal  primal  and  dual  solutions  to  a  perturbed 
problem  of  (P),  where  the  costs  are  perturbed  by  at  most  3y’^/2  and  the  lower  bounds  are  perturbed 
by  at  most  y^^. 


Thus  if  we  choose  e  =  2^0-)  and  y  =  2"®^),  where  L  denotes  the  size  of  the  problem  encoding  in 

binary  (defined  as  in  [1]),  the  homotopy  algorithm  would  terminate  in  steps  with  an 

optimal  primal  dual  solution  pair  to  a  perturbed  problem  of  (P)  and  the  size  of  the  perturbation  is 
2'^).  An  optimal  primal  dual  solution  pair  to  (P)  can  then  be  recovered  by  using,  say,  the 
techniques  described  in  [7]  (also  see  [1]).  Since  the  amount  of  computation  per  step  is  at  most 
0(sv?)  arithmetic  operations  (not  counting  Step  0),  the  homotopy  algorithm  has  a  complexity  of 

0(m3-5.L)  arithmetic  operations.  [We  assume  for  the  moment  that  Step  0  can  be  done  very  "fast". 
See  §5  for  justification.]  This  complexity  can  be  reduced  to  0(m3-3iL)  by  using  Strassen’s 
(impractical)  matrix  inversion  method  [10].  It  may  be  possible  to  reduce  the  complexity  to  ©(m^L) 
by  using  tlie  rank-one  update  technique  described  in  [1],  [11]. 


5.  Algorithm  Initialization 

In  this  section  we  show  that,  for  e  sufficiently  large.  Step  0  of  the  iiomotopy  algorithm  (i.e.  to 
generate  a  primal  dual  pair  (x,u)e  XxSR"  satisfying  p£(x,u)  <  1/2)  can  be  done  very  "fast". 

Suppose  that  (P)  is  in  the  canonical  form  considered  by  Karmarkar  (see  §5  of  [1]  for  details  on 

how  to  transform  general  linear  programs  into  this  canonical  form).  We  claim  that,  for  e  =  2llcll2,  a 

point  (x,u)e  satisfying  p£(x,u)  <  1/2  can  be  found  immediately.  To  see  this,  note  that  in 
Karmarkar's  canonical  form,  A  and  b  have  the  form 
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A' 

0 

A  = 

b  = 

eT 

m 

where  A'  is  some  (n-l)xm  matrix,  and  the  point  e  is  assumed  to  satisfy  Ae  =  b.  Let  x  =  e  and  p  = 
(0 . 0,-1  )T.  Then 

e  +  (ADx)^p  =  e  +  A’*^p 

=  e-e  =  0.  (5.1) 

Hence,  by  (2.7),  x  =  e,  and  the  triangle  inequality, 

P£(x,ep)  =  ll-D^c  +  ee  +  e(ADjj)Tpll2/e 

<  llcll2/e+lle  +  (AD^)Tp|l2 

=  1/2. 


Alternatively,  we  can  solve  the  following  prob.em 

Maximize  Xj  log(xj)  (5.2) 

subject  to  Ax  =  b, 

whose  Karush-Kuhn-Tucker  point  (x,p)  can  be  seen  to  satisfy  (5.1)  (such  a  point  exists  if  the 

feasible  set  for  (P)  is  bounded).  Then,  for  e  =  2110^0112,  we  also  have  Pj(x,ep)  <  1/2.  The  Cunique) 

primal  solution  of  (5.2)  is  sometimes  called  the  analytic  center  [8]  of  the  convex  polyhedron  {  x  I  Ax 
=  b,  X  >  0  ).  Polynomial-time  algorithms  for  solving  (5.2)  are  described  in  [8]  and  [9]  (note  that 
(5.2)  does  not  have  to  be  solved  exactly). 


6.  Conclusion  and  Extensions 


In  this  note  we  have  proposed  a  very  simple  polynomial-time  algorithm  for  linear  programming. 
This  algorithm  solves  a  sequence  of  approximations  to  the  original  problem,  each  augmented  by  a 


logarithmic  penalty  function,  and,  unlike  many  other  Karmarkar-type  algorithms,  uses  no  space 
transformation  of  any  Mnd.  This  algorithm  uses  a  number  of  steps  (namely  0(m*^L))  that  is  the 
lowest  amongst  interior  point  methods  for  linear  programming.  Because  the  primal  and  dual 
solutions  obtained  at  each  step  are  respectively  primal  and  dual  feasible,  determining  termination  is 
particularly  simple:  terminate  if  the  difference  in  their  costs  is  below  a  prespecified  tolerance. 

There  are  many  directions  in  which  our  result  can  be  extended.  Can  the  complexity  be  decreased 
further?  Can  our  algorithm  extend  to  quadratic  programming  or  to  problems  with  upper  bound  on 

the  variables?  Does  there  exist  a  general  class  of  convex  functions  hgtX-^X,  where  X  is  a  convex 

set  in  9?"’,  such  that  hg(x)  — >  ho(x)  pointwise  as  £  i  0  and,  given  an  approximate  minimum  of  hg,  an 

approximate  minimum  of  h„£  (a  is  a  sufficiently  small  scalar  in  (0,1))  can  be  obtained  very 
"quickly"? 
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